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Methods for simulation of biological and/or chemical reac tion pathway, 
biomolecules and nano-molecular systems 

Field of the invention 

The present invention relates to Method for simulation of biological pathways, 
chemical reaction pathways, biomolecules, and nano-molecular-machines by 
using electronic circuits. 

This general field is known as "Molecular Biology" (MB), "Cell Biology" (CB), 
"Nano-technology" (NT), "Chemistry" (CH). When used for analysis of 
disease processes, this field is referred to as "Pathology". When used for the 
examination of the effect of pharmaceutical agents on biological systems, this 
field is referred to as "Drug evaluation" (DE) or "Drug testing" (DT). 

Background of the invention 

The need for simulation of biological pathways and chemical reaction 
pathways by a new approach 

Simulation of the kinetics of biological pathways and chemical reaction 
pathways facilitates the understanding of biological processes, complex 
chemical reactions, disease processes, and the effect of mutations on 
biological systems. It can also be used for testing the effect of 
drugs/molecules on biological systems and disease processes. These are 
described in "Simulation of prokaryotic genetic circuits", H. H. McAddams, and 
A. Arkin. Annu. Rev. Biophys. Biomol. Struct. 27, 199-224 (1998). 

So far, computational methods appear to be the only practically used 
approach for simulation of biological pathways as described in "A 
mathematical model of caspase function in apoptosis", M. Fussenegger, J. E. 
Bailey, and J. Varner. Nature Biotech. 18, 768-774 (2000); "Mathematical 
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modeling of epidermal growth factor receptor signaling through the 
phospholipase C pathway: Mechanistic insights and predictions for molecular 
interventions", J. M. Haugh, A. Wells, D. A. Lauffenburger. Biotech Bioeng 
70, 225-238 (2000); "Computational analysis of biochemical systems, a 
5 practical guide for biochemists and molecular biologists", E. O. Voit. 
Cambridge Press (2000). In these works, the kinetics of binding and reaction 
events in a pathway, under certain approximations, have also been modeled 
by a set of first-order linear differential equations. 

The kinetics of binding and reaction events in a pathway have also been 
described by a set of first-order nonlinear ordinary differential equations as 
proposed in "Systemic properties of ensembles of metabolic networks: 
application of graphical and statistical methods to simple unbranched 
pathways", R. Alves, and M. A. Savageau. Bioinformatics 16, 534-547 (2000); 
and in "Inferring qualitative relations in genetic networks and metabolic 
pathways", T. Akusu, S. Miyano, and S. Kuhara. Bioinformatics 16, 727-734 
(2000). 

While they have been successfully used in studying individual signaling and 
20 metabolic pathways involving dozens of proteins and other molecules, in the 
foreseeable future, computational methods are not expected to have the 
capacity for simulation of large and complex pathways or multiple pathways 
involving hundreds or more number of proteins as described in "Metabolic 
networks: a signal-oriented approach to cellular models" J. W. Lengeler. Biol. 
25 Chem. 381 , 991-920 (2000). 

Cellular events frequently involve "crosstalk" interactions among multiple 
pathways that may include thousands or more proteins and other molecules, 
and the analysis of such complex biochemical networks is further complicated 
30 by the great number of feedback and feedforward loops and of regulatory 
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networks involved in cellular control, as described in "Metabolic networks: a 
signal-oriented approach to cellular models" J. W. Lengeler. Biol. Chem. 381, 
991-920 (2000) and in "Mathematical modeling of epidermal growth factor 
receptor signaling through the phospholipase C pathway: Mechanistic insights 
5 and predictions for molecular interventions", J. M. Haugh, A. Wells, D. A. 
Lauffenburger, Biotech Bioeng 70, 225-238 (2000). 

The kinetics of such a large network is described by up to thousands of first 
order nonlinear differential equations. Significantly more numbers of equations 
10 may be needed for simulation of cellular and multi-cellular systems. 

Given that it is unlikely that the computer resource in the near future is 
capable of solving the "long-time" kinetics of such a large number of nonlinear 
equations, a new approach needs to be introduced. 

15 

The need for structural optimization and molecular dynamics simulation of 
biomolecules and nano-molecular-systems by a new approach 

Molecular dynamics simulation is a widely used method for simulation and 
20 modeling of structures, motions, binding, and thermodynamics of 
biomolecules and nano-molecular-systems. These are described in 
"Computer simulation of biomolecular systems: Theretical and experimental 
applications", Eds. W. F. van Gunsteren, P. K. Weiner, and A. J. Wilkinson. 
ESCOM Leiden (1993); "Computational nanotechnology", R. C. Merkle, 
25 Nanotechnology, 2, 134-141 (1991).; nanostructure engineering as described 
in "Simulated engineering of nanostructures", D.W. Brenner, S.B. Sinnott, 
JA. Harrison, and O.A. Shenderova. Nanotechnology 7 161-167 (1996). It 
has been used in facilitating the study of a variety of important problems such 
as protein folding as described in "Pathways to a protein folding intermediate 
30 observed in a 1 -microsecond simulation in aqueous solution", Y. Duan and 
P.A. Kollman. Science 282, 740-744 (1998); drug design as described in 
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"Binding pathway of retinal to bacterio-opsin: A prediction by molecular 
dynamics simulations", B. Isralewitz, S. Izrailev, and K. Schulten. Biophys. J. 
73, 2972-2979; nanostructure engineering as described in "Simulated 
engineering of nanostructures", D.W. Brenner, S.B. Sinnott, J.A. Harrison, 
5 and O.A. Shenderova. Nanotechnology 7 161-167 (1996). 

Nano-molecular-systems refer to: nano-sized or nano-structured materials 
(such as molecular films, nanotubes, nanoscopic particles with specific 
electronic, optical or magnetic properties), nano-scale molecular mechanical 
10 or manufacturing systems (capable of guiding reactive molecules to 0.1 nm 
precision), other nano-scale devices (molecular motors, carriers, containers, 
pumps, circuits, tools etc.). These are described in "Nanosystems: molecular 
machinery, manufacturing, and computation" by K. Eric Drexler (Wiley 1992). 

15 So far, computational methods are the only practically used approach for 
molecular dynamics simulation as described in "Computer simulation of 
biomolecular systems: Theretical and experimental applications", Eds. W. F. 
van Gunsteren, P. K. Weiner, and A. J. Wilkinson. ESCOM Leiden (1993). 
The molecular dynamics simulation of biomolecules can be described by a set 

20 of second-order nonlinear ordinary differential equations as indicated in "A 
second generation force field for simulation of proteins, nucleic acids, and 
organic molecules", W. D. Cornell et. al., J. Am. Chem. Soc. 117, 5179-5197 
(1995); and in "CHARMM: A program for macromolecular energy, 
minimization, and dynamics calculations", B. R. Brooks, et.al. J. Comp. Chem. 

25 4, 187-217 (1983). While they have been successfully used in studying some 
problems of proteins and other molecules, in the foreseeable future, 
computational methods are not expected to have the capacity for simulation of 
proteins and nucleic acids for orders of magnitude longer than the 
microsecond range. 

30 
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The estimated time scale for folding/unfolding events in proteins, ligand- 
receptor binding/dissociation, and interaction of nano-molecular-machine with 
its substrates ranges from hundreds of milliseconds to tens of seconds as 
described in "Kinetic analysis of the unfolding and refolding of ribonuclease T1 

5 by a stopped-flow double-mixing technique", L. M. Mayr et. al., Biochemistry 
35, 5550-5561 (1996); and in "Studies of the binding of actinomycin and 
related compounds to DNA", W. Muller, and D. M. Crothers, J. Mol. Biol. 35, 
251-290 (1968). In contrast, the longest achieved simulation time is 1 
microsecond for a small protein, consisting of only 36 amino acids, on a 256- 

10 node CRAY T3E supercomputer, as described in "Pathways to a protein 
folding intermediate observed in a 1 -microsecond simulation in aqueous 
solution", Y. Duan and P.A. Kollman. Science 282, 740-744 (1998). Hence 
there is a big gap (5 to 7 orders of magnitude) between simulation time scale 
and that of certain molecular events. Given that it is unlikely the computer 

15 approach in the near future is capable of conducting "long" enough time-scale 
molecular dynamics simulation, a new approach needs to be introduced. 

Summary of the invention 



20 The present invention seeks to address the problems above, and in particular 
to provide a new approach for the simulation of biological and/or chemical 
reaction pathway, and for the molecular dynamic simulation of biomolecules 
and nano-molecular systems by comprising constructing set of binding, 
reaction, concentration and/or molecular dynamics equations and electronic 

25 circuits corresponding to such equations. 

According to a first embodiment, the invention provides a method for 
simulation of at least one biological and/or chemical reaction pathway 
comprising: 
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preparing a map of at least one biological and/or chemical reaction 
pathway; 

constructing at least one set of binding and reaction equation from the 
pathway map; 

5 constructing at least one set of concentration equation for molecules of 

the pathway map; 

constructing an electronic circuit corresponding to every set of 
equation; 

determining simulation of pathway by measuring voltage at two or more 
10 connection points of the circuit. 

In particular, the binding, reaction and/or concentration equation may be a 
linear or non-linear first- or second-order ordinary differential equation (ODE). 

The electronic circuit may comprise at least one of the following circuit units: 
15 linear protein/molecule concentration ±kxj unit; protein/molecule concentration 
power ±kXj a unit; ligand-protein concentration product ±Lxj unit; ligand-protein 
concentration power product ±kLx, a unit; protein-protein concentration product 
±kXiXj unit; protein-protein concentration power product ±kx a Xj b unit; stochastic 
rate constant generator unit that replaces a forward/reverse binding/reaction 
20 rate constant k by a stochastic value; and 

wherein, k is a forward/reverse binding/reaction rate constant, L the 
concentration of a ligand, Xj and x, are the concentration of protein/molecule Xj 
and Xj respectively, a and b are order of power of Xj and Xj. 

25 According to an aspect, the method of the invention further comprises 
maintaining the voltage level of the circuit between two fixed voltage values. 
For example, by: multipling scaling factors to the concentration equation; 
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applying at least one resistor and/or amplifier at one or more connection point 
of the circuit, thereby scaling-down or scaling-up the voltage of one or more 
segment of the circuit; and/or applying automatic gain control circuits. 

5 According to another aspect, the method of the invention further comprises 
adding at least one unit circuit comprising an electronic random noise 
generator and/or a multiplicator amplifier at one or more connection point of 
the circuit. 

10 The method of the invention may further comprise determining the effect of at 
least one drug comprising adding at least one circuit unit associated with a 
receptor protein at one or more connection point of the circuit. 

The method may further comprise determining deficiency, mutation and/or 
15 deletion of at least one protein of the biological pathway, comprising setting a 
voltage ceiling, a voltage range and/or a fixed voltage at one or more 
connection point of the circuit. 

The invention also provides an electronic circuit system for the simulation of at 
20 least one biological and/or chemical reaction pathway, comprising at least one 
electronic circuit representing a set of binding, reaction and/or concentration 
equation. In particular, an electronic circuit system, wherein the binding, 
reaction and/or concentration equation is a linear or non-linear first- or 
second-order ordinary differential equation (ODE). 

25 

The electronic circuit may comprise at least one the following circuit units: 
linear protein/molecule concentration ±kxi unit; protein/molecule concentration 
power ±kxj a unit; ligand-protein concentration product ±l_Xj unit; ligand-protein 
concentration power product ±kLx, a unit; protein-protein concentration product 
30 tkxiXj unit; protein-protein concentration power product ±kXj a Xj b unit; stochastic 
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rate constant generator unit that replaces a forward/reverse binding/reaction 
rate constant k by a stochastic value; and 

wherein, k is a forward/reverse binding/reaction rate constant, L the 
concentration of a ligand, * and Xj are the concentration of protein/molecule x* 
5 and X] respectively, a and b are order of power of x, and X] respectively. 

According to another embodiment, the invention provides a method for 
molecular dynamics simulation of biomolecules and/or nano-molecular 
systems comprising: 

constructing at least one set of equation representing the molecular 
10 dynamics of at least one molecule of the biomolecules and/or the nano- 
molecular systems; 

constructing an electronic circuit representing every set of equation; 

determining molecular dynamics simulation by measuring voltage at 
two or more connection points of the circuit. 

15 

The equation may be a linear or non-linear second order ordinary differential 
equation (ODE). 

According to an aspect, in the method of the invention, the electronic circuit 
20 comprises at least one atom-position circuit unit, wherein the atom-position 
circuit unit represents the position of an atom of a molecule or a molecular 
system. 

In particular, the atom-position circuit unit comprises at least one atom-atom 
interaction circuit subunit, the atom-atom interaction circuit subunit 
25 representing a sub-unit of atom-atom interactions within a molecule or a 
molecular system and comprising internal bond stretch, angle bending, 
torsion, non-bonded units; and/or bond stretch, angle bending, and/or torsion 
units; between at least two nearest sub-units of a molecule. 
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Further, the atom-atom interaction circuit subunit may represent a term in the 
molecular dynamics equation, and wherein the atom-atom interaction circuit 
subunit comprises at least one of the following: bond stretch x unit, bond 
stretch y unit, bond stretch z unit, angle bending x type-A unit, angle bending 

5 x type-B unit, angle bending y type-A unit, angle bending y type-B unit, angle 
bending z type-A unit, angle bending z type-B unit, torsion x type-A unit, 
torsion x type-B unit, torsion y type-A unit, torsion y type-B unit, torsion z type- 
A unit, torsion z type-B unit, non-bonded x unit, non-bonded y unit, non- 
bonded z unit, hydrogen-bond x unit, hydrogen-bond y unit, and hydrogen- 

10 bond z unit; and wherein x, y, and z represent the coordinates of each atom of 
the molecule, and type-A represents the case of the atom being in the middle- 
position of an angle bending or torsion connection with other atoms, and type- 
B represents the case of the atom being in the end-position of an angle 
bending or torsion connection with other atoms. 

15 The method may further comprise maintaining the voltage level in the circuit 
between two fixed voltage values. For example, by wherein x, y and z 
represent the coordinate of the molecule, and the voltage level of the circuit is 
maintained between two fixed voltage values by: 

applying scaling factors to the x, y and z coordinates and to the 
20 molecular dynamic equation; 

applying at least one resistor and/or amplifier at one or more 
connection point of the circuit, thereby scaling-down or -up the voltage of one 
or more segment of the circuit; and/or 

applying automatic gain control circuits. 

25 The invention also provides a circuit group representing the interaction pattern 
in the chemical structure of a molecule or a sub-unit of interaction pattern in 
the chemical structure of a molecule comprising: 
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a bond stretch connection between each atom pair of the molecule 

covalently bonded to each other; 

an angle bending connection pair between a first atom and other two 

atoms; 

5 a torsion connection bundle between a first atom and other three 

atoms; and 

a non-bonded connection between each atom pair whose atoms are at 
least four bonds away from each other. 

10 The electronic circuit may comprise at least one circuit unit, the circuit unit 
comprising at least one circuit group, the circuit group representing a sub-unit 
of interaction pattern in the chemical structure of a molecule and comprising 
internal bond stretch, angle bending, torsion, non-bonded units; and/or bond 
stretch, angle bending, and/or and torsion units; between at least two nearest 

15 sub-unit of a molecule. 

The circuit unit may represent a term in the molecular dynamic equation, and 
wherein the circuit unit comprises at least one of the following: bond stretch x 
unit, bond stretch y unit, bond stretch z unit, angle bending x type-A unit, 
20 angle bending x type-B unit, angle bending y type-A unit, angle bending y 
type-B unit, angle bending z type-A unit, angle bending z type-B unit, torsion x 
type-A unit, torsion x type-B unit, torsion y type-A unit, torsion y type-B unit, 
torsion z type-A unit, torsion z type-B unit, non-bonded x unit, non-bonded y 
unit, non-bonded z unit, hydrogen-bond x unit, hydrogen-bond y unit, and 
25 hydrogen-bond z unit; and wherein x, y, and z represent the coordinates of 
each atom of the molecule, and type-A unit represents the case of the atom 
being in the middle-position of an angle bending or torsion connection with 
other atoms, and type-B represents the case of the atom being in the end- 
position of an angle bending or torsion connection with other atoms. 
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The invention further relates to a method for the manufacture of an electron* 
circuit representing at least one biomolecule and/or nano-molecular system, 
the electronic circuit comprising at least a unit circuit comprising at least a 
circuit group, wherein the circuit group represents a sub-unit of interacts 
5 pattern in the chemical structure of a molecule or a molecular system 
comprising: 

introducing a bond stretch between each atom of a pair of atoms 

covalently bonded to each other; 

introducing an angle bending connecting pair between a first atom and 

10 other two atoms; 

introducing a torsion connection bundle between a first atom and other 

three atoms; and 

introducing a non-bonded connection between each atom pair whose 
atoms are at least 4 bonds away from each other. 
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Brief description of th e drawings 

Figure 1 is a flow diagram that illustrates one embodiment of a process for 
simulating a single or multiple biological or chemical reaction pathways by 
electronic circuits (100). 

Figure 2 is the map of an illustrative biological pathway, part of the caspase 
pathway. 

Figure 3A, 3B are the maps of the actual circuits for an illustrative biological pathway, 
part of the caspase pathway, wherein Figure 3A represents the circuit for Equation I 
{!b. = klXi - k> .L Xi+ o> l ) an d F ® ure 38 represents the circuit for Equation II 

dt 

(^ = fc 1+ Lx, -(*,_+ fc 2 > 2 )- 
at 
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Figure 4 is a map representing an electronic emulator for the entire pathway 
governed by twenty ODEs, of the caspase pathway, where blocks Equation 1 
to 20 represent the circuits that emulate the twenty ODEs. 

5 Figure 5 is the measured voltage changes, representing variations of protein 
concentrations in the illustrative biological pathway, from the designed circuit 
for this pathway. 

Figure 6 is a schematic diagram of the construction of a section of electronic 
10 circuits representing an ODE associated with a binding/reaction event in a 
biological or chemical pathway. Note that switches 11 and 12 are closed at t=0 
when the initial condition is Xi=a or dxi/dt=B. Each switch can also be used 
when Xi or dxj/dt is set at a fixed value or a fixed range. 

15 Figure 7 is a schematic diagram of the circuit for the first ODE of the 
illustrative caspase pathway. Note that switches 11 and 12 are closed at t=0 
when the initial condition is Xi=a or dXj/dt=B. Each switch can also be used 
when x or dx/dt is set at a fixed value or a fixed range. 

20 Figure 8 is a schematic diagram of the design plans for the special circuit 
units representing each term in the ODE of a binding/reaction event in a 
pathway. 

Figure 9 is a schematic diagram of the design plans for the special circuit 
25 units representing each term in the ODE of a binding/reaction event in a 
pathway. The random number £, is in the range: 0<^<1. 
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Figure 10 is a schematic diagram of the construction of a section of electronic 
circuits representing a molecular dynamics simulation ODE for the x- 
coordinate of an atom. 
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Figure 1 1 is a schematic diagram of the construction of a section of electronic 
circuits representing a molecular dynamics simulation ODE for the y- 
coordinate of an atom. 

Figure 12 is a schematic diagram of the construction of a section of electronic 
circuits representing a molecular dynamics simulation ODE for the z- 
coordinate of an atom. 

Figure 1 3 is a schematic diagram of the bond stretch connections of the circuit 
group representing a subunit, amino acid, of an illustrative biomolecule, 
protein. 

Figure 14 is a schematic diagram of the angle bending connections of the 
15 circuit group representing a subunit, amino acid, of an illustrative biomolecule, 
protein. 

Figure 15 is a schematic diagram of the torsion connections of the circuit 
group representing a subunit, amino acid, of an illustrative biomolecule, 
20 protein. 

Figure 16 is the structure of an illustrative nano-molecular-machine, a 
rudimentary mimic of the enzyme ribonuclease, which has been constructed 
from the bucket-shaped cyclodextrin molecule (a), midazole functional groups 
25 on the rim (X in (b)) help to catalyze the region specific hydrolysis of a cyclic 
phosphodiester (b), which is bound in the hydrophobic cavity. 

Figure 17 is a schematic diagram of the bond stretch connections of the circuit 
group representing a subunit of an illustrative nano-molecular-machine, a 
30 rudimentary mimic of the enzyme ribonuclease. 
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Figure 18 is a schematic diagram of the angle-bending connections of the 
circuit group representing a subunit of an illustrative nano-molecular-machine, 
a rudimentary mimic of the enzyme ribonuclease. 

5 

Figure 19 is a schematic diagram of the torsion connections of the circuit 
group representing a subunit of an illustrative nano-molecular-machine, a 
rudimentary mimic of the enzyme ribonuclease. 

10 Figure 20 is a schematic diagram of the bond stretch x unit, representing the 
x-component of a bond stretch term. 

Figure 21 is a schematic diagram of the bond stretch y unit, representing the 
y-component of a bond stretch term. 

15 

Figure 22 is a schematic diagram of the bond stretch z unit, representing the 
z-component of a bond stretch term. 

Figure 23 is a schematic diagram of the angle bending x type-A unit, 
20 representing the x-component of a bond angle bending term in which the 
atom is in the middle. 

Figure 24 is a schematic diagram of the angle bending x type-B unit, 
representing the x-component of a bond angle bending term in which the 
25 atom is at one end. 

Figure 25 is a schematic diagram of the angle bending y type-A unit, 
representing the y-component of a bond angle bending term in which the 
atom is in the middle. 
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Figure 26 is a schematic diagram of the angle bending y type-B unit, 
representing the y-component of a bond angle bending term in which the 
atom is at one end. 

Figure 27 is a schematic diagram of the angle bending z type-A unit, 
representing the z-component of a bond angle bending term in which the 
atom is in the middle. 

Figure 28 is a schematic diagram of the angle bending z type-B unit, 
representing the y-component of a bond angle bending term in which the 
atom is at one end. 

Figure 29 is a schematic diagram of the torsion x type-A unit, representing the 
x-component of a bond torsion term in which the atom is in the middle. 

Figure 30 is a schematic diagram of the torsion x type-B unit, representing the 
x-component of a bond torsion term in which the atom is at one end. 

Figure 31 is a schematic diagram of the torsion y type-A unit, representing the 
20 y-component of a bond torsion term in which the atom is in the middle. 

Figure 32 is a schematic diagram of the torsion y type-B unit, representing the 
y-component of a bond torsion term in which the atom is at one end. 

25 Figure 33 is a schematic diagram of the torsion z type-A unit, representing the 
z-component of a bond torsion term in which the atom is in the middle. 

Figure 34 is a schematic diagram of the torsion z type-B unit, representing the 
z-component of a bond torsion term in which the atom is at one end. 
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Figure 35 is a schematic diagram of the non-bonded x unit, representing the 
x-component of a non-bonded interaction term. 

Figure 36 is a schematic diagram of the non-bonded y unit, representing the 
5 y-component of a non-bonded interaction term. 

Figure 37 is a schematic diagram of the non-bonded z unit, representing the 
z-component of a non-bonded interaction term. 

1 o Figure 38 is a schematic diagram of the H-bond x unit, representing the x- 
component of a hydrogen-bond term. 

Figure 39 is a schematic diagram of the H-bond y unit, representing the y- 
component of a hydrogen-bond term. 

15 

Figure 40 is a schematic diagram of the H-bond z unit, representing the z- 
component of a hydrogen-bond term. 

Figure 41 is a schematic diagram of the rABC circuit unit and the AqA circuit 
20 unit. 

Figure 42 is a schematic diagram of the Ua circuit unit, ss circuit unit, and si 
circuit unit. 

25 Figure 43 is a schematic diagram of the pq circuit unit. 



Figure 44 is a schematic diagram of the <j> circuit unit and the sin<j> circuit unit. 
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Detailed description o f the invention 



Method for simulation of biological and/or chemical reaction pathway 

5 According to a first embodiment, the present invention provides a method for 
simulation of a single or multiple biological and/or chemical reaction' 
pathways. 

It has particular application in simulation of large networks of pathways that 
10 are beyond the scope of computer methods. 

It provides a potentially useful means for facilitating the study of cellular and 
multi-cellular processes, disease processes, effect of gene or protein 
mutations on biological systems, and effect of drugs on biological systems. 

15 

A single or a set of biological pathways and/or chemical reaction pathways 
can be represented by a map composed of a network of boxes or labels 
(reaction objects) connected by arrows (binding/reaction paths). This kind of 
representation is widely used in popular biological pathway databases such 
20 as SPAD (htt p://www.q rt kyj iRhii-u.ac.ip/spad/) and KEGG 
mttn-//www.aen" " 0 ad j P /kenn/keaa2.htmlV and as described in 
"Mathematical modeling of epidermal growth factor receptor signaling through 
the phospholipase C pathway: Mechanistic insights and predictions for 
molecular interventions", J. M. Haugh, A. Wells, D. A. Lauffenburger. Biotech 
25 Bioeng 70, 225-238 (2000). Each binding/reaction event in these pathways 
can be described by a binding/reaction equation, whose form depends on 
whether the event is reversible or irreversible. 

The set of binding and/or reaction (binding/reaction) equations for pathways 
30 can be modeled by a set of linear and/or non-linear first-order and/or second- 
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order ordinary differential equations. For example, non-linear first-order 
ordinary differential equations (ODEs) derived from the balance equation and 
generalized mass-action kinetics as described in "Understanding the control 
of metabolism" D. Fell, Portland Press, London, UK (1996) and in "The 

5 regulation of cellular systems" R. Heinrich, and S Schuster, Chapman & Hall, 
New York (1996). Each ODE describes the time-dependent behaviour of the 
concentration of a protein/molecule in the pathways. Further, equations 
representing the concentration of the molecules of the pathway map of the 
biological and/or chemical reaction pathway are also provided. The 

10 concentration equations can also be at least a set of linear and/or non-linear 
first-order and/or second-order ordinary differential equations. For example, 
non-linear first-order ordinary differential equations (ODEs). 

According to this embodiment of the invention, at least one set of electronic 
15 circuits represents the at least one set of the above equations. In particular, 
binding, reaction and/or concentration non-linear ODEs. 

Electronic circuits have been proposed to solve certain first-order and second- 
order nonlinear ordinary differential equation as described in "The transition 
20 from solitons to chaos in the solution of the logistic equation" M. I. Sobhy, and 
A. S. Burman. Int. J. Bifurcation Chaos 10, 2823-2829 (2000) as well as linear 
ordinary differential equation (up to second-order) as described in "The analog 
computer as an aid to the teaching elementary quantum mechanics", M. K. 
Summers. Phys. Educ 13, 22-27 (1996). 
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However, electronic circuits have never been disclosed nor suggested for the 
purpose of simulation of biological and/or chemical reaction pathway. 



Operational amplifiers, function circuits and other electronic components can 
30 be used to construct operations of differentiation, integration, sum, 
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subtraction, multiplication, inversion, exponential, logarithm, power, and 
others as described in "Electronic circuits and design" D. A. Neamen. McGraw 
Hill (2001), "Function circuits: Design and applications" Y.J. Wong and W.E. 
Ott, McGraw-Hill (1976), "Modern operational circuit design" J.I. Smith, Wiley- 

5 Interseicne (1971). Hence, a set of nonlinear or linear ordinary differential 
equations composed of these terms can be represented by a set of electronic 
circuits composed of these operational amplifiers and other electronic 
components. The time-dependent behaviour of the concentration of a 
protein/molecule can be measured by the voltage at a specific point in the 

10 circuits. 

Accordingly, for simulation of biological and/or chemical reaction pathways, 
the time-dependent behavior of the concentration of one or more reaction 
objects, defined as proteins/nucleic 

15 acids/ligands/substrates/reactants/products, in the pathways can be derived 
by measuring the voltage at specific points in the circuits. 

In order to accurately describe the microscopic kinetics of binding or reaction 
events in pathways, all possible states for each protein or molecule need to 

20 be considered. Examples of these states include active, inactive, 
ligand/substrate-bound, dimmer, and any other state. The corresponding 
circuit section for the concentration in a particular state is typically composed 
of up to a few dozen operational amplifiers, a similar number of connection 
points (two for measuring its voltage, a few for connection to external voltage 

25 sources, and several for connection to and from other circuit units), and a 
number of small components such as switches etc. The circuits for a large 
pathway or multiple pathways may thus contain up to thousands or more such 
sections connected to each other and to the external voltage sources, which 
is feasible to construct. 
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Accordingly, the first embodiment provides a method for simulation of at least 
one biological and/or chemical reaction pathway comprising: 

preparing a map of at least one biological and/or chemical reaction 
pathway; 

5 constructing at least one set of binding and reaction equation from the 

pathway map; 

constructing at least one set of concentration equation for molecules of 
the pathway map; 

constructing an electronic circuit corresponding to every set of 
10 equation; 

determining simulation of pathway by measuring voltage at two or more 
connection points of the circuit. 

The binding, reaction and/or concentration equation may be a linear or non- 
15 linear first- or second-order ordinary differential equation (ODE). Preferably, 
the binding, reaction and/or concentration equation is a non-linear first-order 
ordinary differential equation (ODE). 

In a further aspect of this invention, methods are disclosed for maintaining the 
20 voltages in the circuits (concentrations of molecules and functions of 
concentrations) within the allowed range. Scaling factors (which are non-zero 
real numbers) are applied to the concentrations and to the kinetic equations 
so as to keep these concentrations and each term in the equations within the 
allowed range. In addition, resistors and amplifiers can be used to scale-down 
25 and scale-up the voltages at certain segment of the circuits. Moreover 
automatic gain control circuits can also be used to ensure the voltages are 
within the required range. The automatic gain control circuits are described in 
"Function circuits: Design and applications" Y.J. Wong and W.E. Ott, 
McGram-Hill, 1976. 
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Accordingly, the invention provides a method further comprising maintaining 
the voltage level of the circuit between two fixed voltage values. 

5 In particular, the voltage level of the circuit is maintained between two fixed 
voltage values by: 

multipling scaling factors to the concentration equation; 
applying at least one resistor and/or amplifier at one or more 
connection point of the circuit, thereby scaling-down or scaling-up the voltage 
10 of one or more segment of the circuit; and/or 
applying automatic gain control circuits. 

In a further aspect of this invention, stochastic effect on the binding/reaction 
events, which is important in cases involving low concentration or short-time 

15 event, can be taken into account by replacing the relevant binding/reaction 
rate constants by stochastic values (random variations around average rate 
constants determined by experiments or by theoretical methods). This can be 
accomplished in the circuits by adding at specific sites special circuit units 
composed of an electronic random noise generator device, such as Wayne's 

20 random noise generators, and a multiplication amplifier. 

Accordingly, the method of the invention further comprises adding at least one 
unit circuit comprising an electronic random noise generator and/or a 
multiplicator amplifier at one or more connection point of the circuit. 

25 

In a further aspect of the invention, design plans for several circuit units are 
introduced to represent one or more, or even all possible terms in a circuit 
equation. Accordingly, in the method of the invention, the electronic circuit 
may comprise at least one of the following circuit units: linear protein/molecule 
30 concentration ±kx, unit; protein/molecule concentration power ±kXi a unit; 
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ligand-protein concentration product ±LXi unit; ligand-protein concentration 
power product ±kLXj a unit; protein-protein concentration product ±kXiXj unit; 
protein-protein concentration power product ±kXi a Xj b unit; stochastic rate 
constant generator unit that replaces a forward/reverse binding/reaction rate 
5 constant k by a stochastic value; and 

wherein, k is a forward/reverse binding/reaction rate constant, L the 
concentration of a ligand, * and Xj are the concentration of protein/molecule * 
and Xj respectively, a and b are order of power of * and Xj. 

10 By introducing additional circuit sections or by modifying the circuits, it is 
possible to simulate events such as the effect of drugs, protein or gene 
mutations, and diseases on a pathway or a group of pathways. 

Therefore, in a further aspect of this invention, effect of drugs can be 
15 simulated by adding one or more circuit units (the design plan for each of 
these units is disclosed) at specific sites in the circuits that are associated with 
the respective receptor proteins. These added units represent the forward and 
reverse binding/reaction process between these drugs and their respective 
receptor proteins. 

20 

Accordingly, the method of the invention, further comprises determining the 
effect of at least one drug comprising adding at least one circuit unit 
associated with a receptor protein at one or more connection point of the 
circuit. 

25 

In a further aspect of this invention, deficiency or mutation of one or more 
proteins can be simulated by setting a voltage ceiling, or a voltage range, or a 
fixed voltage value at specific points in the circuits. 
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Accordingly, the method of the invention further comprises determining 
deficiency, mutation and/or deletion of at least one protein of the biological 
pathway, comprising setting a voltage ceiling, a voltage range and/or a fixed 
voltage at one or more connection point of the circuit. 

5 

The method of claim 1 , wherein the biological pathway comprises at least one 
object, and wherein the object is protein, nucleic acid, ligand, substrate, 
inhibitor or antagonist, activator or agonist, reactant and/or reaction product. 

10 According to the first embodiment of the invention, it is also provided an 
electronic circuit system for the simulation of at least one biological and/or 
chemical reaction pathway, comprising at least one electronic circuit 
representing a set of binding, reaction and/or concentration equation. 

15 In the electronic circuit system, the binding, reaction and/or concentration 
equation may be a linear or non-linear first- or second-order ordinary 
differential equation (ODE). Preferably, a non-linear first-order ODE. 

According to a particular aspect, the invention relates to the disclosure of 
20 design plans for the construction of circuit units that represent one or more or 
even all possible terms in the ODEs describing a biological or a chemical 
reaction pathway. 

Accordingly, in the electronic circuit system of the invention, the electronic 
25 circuit may comprise at least one, or even all, the following circuit units: linear 
protein/molecule concentration ±kx, unit; protein/molecule concentration power 
±kXi a unit; ligand-protein concentration product ±LXi unit; ligand-protein 
concentration power product ±kl_Xi a unit; protein-protein concentration product 
±kXjXj unit; protein-protein concentration power product ±kXj a Xj b unit; stochastic 
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rate constant generator unit that replaces a forward/reverse binding/reaction 
rate constant k by a stochastic value; and 

wherein, k is a forward/reverse binding/reaction rate constant, L the 
concentration of a ligand, x* and Xj are the concentration of protein/molecule Xj 
and Xj respectively, a and b are order of power of Xj and Xj respectively. 

Referring now to the drawings and, in particular, to Figure 1, there is disclosed 
a process 100 for carrying out a method for the simulation of a single or 
multiple biological or chemical reaction pathways. In the process 100, 
simulation of the pathways is conducted on a specifically designed set of 
electronic circuits by the procedure as described below and given in the flow 
chart of Figure 1. 

The process 100 begins at a state 110, wherein the pathways are presented 
as a map of a network of boxes or labels (proteins or molecules) connected 
by arrows (binding/reaction paths). Figure 2 shows the map of a part of 
caspase pathway as an illustration. 

The process 100 then moves to a sate 120 wherein the binding or reaction 
equation for every process in the pathway is constructed from the pathway 
map. Based on whether it is reversible or irreversible, the equation of an 
individual binding /reaction step takes the form: 

Reversible binding/reaction: Ri + R2 + . . . \— v Pi + P2 + . . . 

K- 

K* 

Irreversible binding/reaction: Ri + R2 + . . . x Pi + P2 + . . . 

where Ri, R2, ... are binding/reaction objects (proteins/molecules etc) and P 1f 
P 2 , ... are product objects (proteins/molecules etc), K+ and K. are forward and 
reverse binding/reaction rate constant respectively. 
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20 



The first few reaction equations for the illustrative caspase pathway in Figure 
2 are given as: 



K3» 



Event 1: Xi + L v= v Xz Event 3: Xs + X4 v= x X 5 



ki- 



JO- 



10 Event 2: X 2 ---- x X 3 Event 4: X4 + Xs v== v Xe 



K4- 



where l = Fas ligand, Xi=ligand-free Fas surface receptor, x 2 =ligand-bound 
Fas surface receptor, x 3 =clustered ligand-bound Fas surface receptor, 

15 X4=FADD protein, x 5 =complex of FADD-Fas receptor and k i+ and k,. are 

the forward and reverse reaction rate constant in the i-th path respectively. 
These rate constants can either be obtained from experimental data, or from 
theoretical computations, or estimated by empirical means based on statistical 
analysis of available experimental data of available of relevant reactions. 



In general, the reaction equation for the i-th path is: 



25 Path i: S j a u Xj \=== 2 k b ik x k 

k,. 



where a q and b ik are the stoichiometric coefficient of the j-th and k-th reaction 

30 objects respectively; j=1 N,; k=1, ...,N k ; and Nj and N k are number of 

reaction objects and number of product objects respectively. 

The process 100 then moves to a state 130 where in the corresponding 
kinetic equations for the concentration of each protein or molecule are derived 
35 from the balance equation and generalized mass-action kinetics as described 
in "Understanding the control of metabolism" D. Fell, Portland Press, London, 
UK (1996) and in 'The regulation of cellular systems" R. Heinrich, and S 
Schuster, Chapman & Hall, New York (1996). Each equation is a non-linear 
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first order ordinary differential equation (ODE) of the concentration of a 
protein/molecule in the pathways. In general, for a set of equations involving 
the n-th reaction object: 

5 Ki* 

a,„ Xn + 2 j a tj Xj \— v S k b lk x k 

Kl 

Ki + 

2 j- a ir Xj' \== v ar„ Xn + 2 k - b rk - x k> 
10 Kr- 

the corresponding ODE is: 
1 5 dx„/dt = k,. n k (x k ) b,k + kr. n j. (x r ) air + . . . 

. k, + (x„) ain n , (Xj) aij - kr. (x n ) aln n k . (x k .) blk ' -... + «>„ 

where x n is the concentration of object x n , co n is its rate of 
20 expression/production, ay and by are the power of the concentration in cases 
where the reaction rate is determined by power law representation. 

The first few ODEs for the illustrative caspase pathway are: 

25 

dx^dt = ki. Xa - k 1+ L + oo, dxj/dt = k 2 Xz + k 3 . Xs - k 3+ x 3 x 4 

dx 2 /dt= k 1+ LXi - (k,.+ k 2 )x 2 dxjdt= k 3 . Xs + k 5 . Xe - k 3 * x 3 X4- k s+ X4X S +<d 4 

where the ligand concentration L and rate of expression a> n of protein x n are 
30 assumed to be time-dependent function. 

Referring back to Figure 1, the process 100 then moves to a state 140 
wherein a network of electronic circuits are constructed to represent the ODEs 
of a pathway or a network of pathways. This circuit network is composed of 
35 specially designed circuit units as disclosed in this patent and known 
components such as operational amplifiers, voltage sources representing 
time-dependent distribution of ligands or molecules, and other electronic 
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circuit components such as switches and voltage sources. Each equation has 
a circuit, which is connected to other circuits. An equation normally contains 
no more than a dozen terms or variables, and thus its circuit needs to have no 
more than a dozen connection points to other circuits. Each circuit starts with 
5 a point x,, which is connected to a switch that controls the initial condition, or 
the range, or the fixed value of voltage. This point is connected to other parts 
of the circuit by a combination of one or more of the disclosed specially 
designed circuit units together with a differentiator and a summing amplifier as 
illustrated in Figure 6. Figure 7 shows the actual circuit design for the first 
10 ODE of the illustrative caspase pathway. The differentiator, summing 
amplifier, and inverting amplifier are described in "Electronic circuit analysis 
and design, 2 nd edition" D.A. Neamen, McGraw-Hill, 2001. The multiplier is 
described in "Fundamentals of linear circuits" T.L. Floyd, Macmillan Publishing 
Co. 1992. 



15 



20 



25 



The process 100 then moves to a simulation state 150 wherein the time- 
dependent kinetics of the concentration of one or more proteins or molecules 
in a pathway or a group of pathways can be determined by measurement of 
voltages at specific sites of the circuits. For the circuit illustrated in Figure 6, 
the voltage at x ( corresponds to the concentration of protein/molecule x,. For 
the circuit illustrated in Figure 7, the voltage at Xi corresponds to the 
concentration of protein xi. 

Construction Plan for the Specially Designed Circuit Units Representing 
Different Terms in an ODE For Pathway Simulation 



Referring now to Figure 8 and Figure 9, there are disclosed design plans for 
30 construction of the special circuit units according to the invention representing 
all possible terms of an ordinary differential equation associated with a 
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biological and/or chemical reaction pathway. The operational amplifiers and 
random noise generator in these units are all described in the literature. The 
inverting amplifier and noninverting amplifier are described in "Electronic 
circuit analysis and design, 2 nd edition" D.A. Neamen, McGraw-Hill, 2001. The 
multiplier is described in "Fundamentals of linear circuits" T.L. Floyd, 
Macmillan Publishing Co. 1992. The exponentiator is described in "Function 
circuits: Design and applications" Y.J. Wong and W.E. Ott, McGram-Hill, 
1976. The random noise generator is described in "A Noise-Based IC 
Random Number Generator for Applications in Cryptography," Petrie, C.S. 
and Connelly, J.A., IEEE Transactions on Circuits and Systems, Part I, Vol. 
47, May 2000, pp. 615-621. 

Methods for Maintaining the Voltages in the Circuits Within Allowed Range 

Normally the concentration of a molecule can be scaled into the allowed 
voltage range of a circuit. In case that this may not be so, or in case that the 
derivative of concentration exceeds the voltage range, or in case that some 
terms in the kinetic equation may exceed the voltage range, the following 
methods can be used to maintain the voltage level of the circuit between two 
voltage values: 

(1) Applying scaling factors to the concentrations or to the pathway 
simulation ODEs so as to keep the concentrations and each term in 
the equations within the allowed range. 

This can be accomplished by multiplying an appropriately selected 
small number 5 to the pathway simulation ODEs (that is, using non- 
zero real numbers) 
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(2) Resistors are used to scale-down the voltages at the input end of a 
segment of the circuits and amplifiers are used at the output end of the 
segment to sale-up the voltages to the proper level. 

This can be accomplished by applying at least one resistor and/or 
amplifier at one or more connection point of the circuit, thereby scaling- 
down or scaling-up the voltage of one or more segment of the circuit. 

(3) Automatic gain control circuits can be added to the circuits to ensure 
the voltages are within the required range. 

Simulation of Effect of Drugs 

The effect of drugs on a pathway or a group of pathways can be simulated by 
adding ligand-protein concentration product units to the circuit sections of the 
drug receptor proteins. This corresponds to case 3 in Figure 6. The ligand- 
protein concentration product units are shown in Figure 9. 

Illustrative example: Electronic pathway emulator of a caspase pathway 

As mentioned before in the description of the first embodiment of the 
invention, a pathway can be described by a set of binding/reaction equations. 
This set of binding and/or reaction equations can be further modeled by a set 
of ODEs. An illustrative biological pathway is part of caspase pathway in 
Figure 2 is governed by the following 20 simultaneous ODEs. 
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(0 *i = * x 1 -k lt Lx l +<o i 

dt 1 

(II) ^ = + *,)*, 
5 * 

(III) ^z^+^-MA 

a/ 

(IV) = k z _x s + * 5 _x 6 - * 3+ * 3 x 4 - * 5+ * 4 *s + G> 4 
dt 

10 (V) ^ = k s-*6 + ^ Xi x 4 - k,_x 5 - k s+ x,x 5 

(vii) = k <- x * + * 8 -* 9 " * tf ** s * 7 ~ * 8+ * 7 * 8 + 6,7 
15 (viii) ^ = * 8+ * 9 + kx l ~ Xn + K * x&Xl " * 6 '* 8 " * 8 ~* 7 * 8 " * ll+ * 8 * 10 



(IX) 

(X) ^ = + fc H-*U - *1**«*I0 - *lf*.*H> + ®10 



20 (XI) ^u- = jfc 10+ ^ 10 -V^ 

(XI I ) 12 = ^n+^^io ~ ^ll-'u 

a/ 

(X|||) = fc 9 JC 9 + *15-*1S + *15*15 + *17-*M - *i3 + *13*14 " *17 + *13*17 

dx 

(xv) iir = /Ci3 ^ i3Xi4 ~ ( * 13 - + * i5)Xis 
25 (xvi ) = fc,5 * ls + *"- JC20 " kl ° M ' 

(XVII) ^-^-*^ + -* 

^13 _ K y v —It X 
— — — ^17+^13^17 "-17-^18 

(XVIII) * 

(xix) f 

— ^19+^16^19 ^19-^20 

(XX) * 

30 
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5 where the ligand concentration L and rate of expression co n can change with 
time (a time-dependent function). Variable x1 to x20 represent the protein 
concentration at the respective paths. This set of nonlinear ODEs can be 
represented by a corresponding set of electronic circuits composed of 
different function blocks. Thus, an electronic biological pathways emulator 

10 can be built to mimic and study the biological process. 

Among the given 20 simultaneous ODEs, Eqn. (I) to (IV) are of typical types. 
These typical equations can be emulated using the circuits shown in Figures 
3A, 3B and in Figure 4. 

15 

The emulator circuit has been constructed and the measured voltages, 
representing variations of protein concentrations in the illustrative biological 
pathway, are shown in Figure 5, where voltage x1, x2, ... x20 represent the 
concentration of each protein in a particular state (active or inactive). The 
20 results have shown the expected behaviour as found in computer simulation. 

Molecular dynamics simulations of biomolecules and/or nano-molecular- 
system 

25 

According to a second embodiment, the invention relates to molecular 
dynamics simulations of biomolecules and/or nano-molecular-systems. 

Nano-molecular-systems refer to: nano-sized or nano-structured materials 
30 (such as molecular films, nanotubes, nanoscopic particles with specific 



32 



electronic, optical or magnetic properties), nano-scale molecular mechanical 
or manufacturing systems (capable of guiding reactive molecules to 0.1 nm 
precision), other nano-scale devices (molecular motors, carriers, containers, 
pumps, circuits, tools etc.). These are described in "Nanosystems: molecular 
5 machinery, manufacturing, and computation" by K. Eric Drexler (Wiley 1 992). 

The molecular dynamics of a protein/nucleic acid/organic molecule can be 
represented by a set of linear and/or non-linear second-order ordinary 
differential equations (ODEs) derived from the Newton's second law of motion 

10 given the energy functions, force fields and a starting 3D structure of the 
system, as described in "A second generation force field for simulation of 
proteins, nucleic acids, and organic molecules", W. D. Cornell et. al., J. Am. 
Chem. Soc. 117, 5179-5197 (1995); and in "CHARMM: A program for 
macromolecular energy, minimization, and dynamics calculations", B. R. 

15 Brooks, et.al. J. Comp. Chem. 4, 187-217 (1983). Each ODE describes the 
time-dependent change of the x- or y- or z- coordinate of each atom in a 
studied molecule. 

According to the second embodiment of the invention, a corresponding set of 
20 electronic circuits is proposed to represent this set of linear or non-linear 
ODEs. This set of electronic circuits is composed of operational amplifiers, 
function circuits and other electronic circuit elements. A combination of one or 
more specially designed circuit units (composed of operational amplifiers, 
function circuits and other electronic components together with differentiator 
25 and summing amplifier) is used to represent each term of the ODEs. The 
time-dependent position of each atom in a protein/nucleic acid/organic 
molecule can be measured by the voltage at a specific point in the circuits. 

Electronic circuits have been proposed to solve certain first-order and second- 
30 order nonlinear ordinary differential equation as described in "The transition 
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from solitons to chaos in the solution of the logistic equation" M. I. Sobhy, and 
A. S. Burman. Int. J. Bifurcation Chaos 10, 2823-2829 (2000) as well as linear 
ordinary differential equation (up to second-order) as described in "The analog 
computer as an aid to the teaching elementary quantum mechanics", M. K. 
Summers. Phys. Educ 13, 22-27 (1996). However, electronic circuits have 
never been disclosed nor suggested for the purpose of molecular dynamics 
simulations of biomolecules and/or nano-molecular systems. 

If a 3D structure is not available, a starting structure of a molecule can be 
generated from its molecular bond connectivity profile and from the standard 
average bond length, angle and torsion. Solution of these molecular dynamics 
simulation equations gives the time evolution of the position of each atom of a 
studied molecule. All the terms in this equation can be represented by a 
combination of existing operational amplifiers, function circuits and other 
electronic components. Hence it is feasible to use electronic circuits to solve 
these equations and thus conduct molecular dynamics simulations. 

A protein normally consists of thousands of atoms. Each atom has three 
equations describing its x-, y-, or z- coordinate respectively. The 
corresponding circuit section for each coordinate is typically composed of up 
to a few dozen bond-stretch/angle/torsion units plus up to thousands of non- 
bonded units. If a force field with separate hydrogen bond terms is used, then 
this circuit section also contains a few additional units for hydrogen bonds. If 
explicit solvent is used, then this circuit section also contains additional non- 
bonded units related to the non-bonded interactions with water molecules. 
Each unit consists of approximately a few dozen operational amplifiers and 
function amplifiers, a similar number of connection points (two for measuring 
its voltage, a few for connection to external voltage sources, and the rest for 
connection to and from other circuit units), and a number of small components 
such as switches etc. The circuits for a typical protein may thus contain up to 
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thousands or more such sections connected to each other and to the external 
voltage sources, which is feasible to construct. 

Operational amplifiers, function circuits and other electronic components can 
5 be used to construct operations of differentiation, integration, sum, 
subtraction, multiplication, inversion, exponential, logarithm, power, and 
others as described in "Electronic circuits and design" D. A. Neamen. McGraw 
Hill (2001), "Function circuits: Design and applications" Y.J. Wong and W.E. 
Ott, McGraw-Hill (1976), "Modern operational circuit design" J.I. Smith, Wiley- 
10 Interseicne (1971). Hence, a set of nonlinear or linear ordinary differential 
equations composed of these terms can be represented by a set of electronic 
circuits composed of these operational amplifiers and other electronic 
components. The time-dependent behavior of the concentration of a 
protein/molecule can be measured by the voltage at a specific point in the 
15 circuits. 

The present invention provides a method for molecular dynamics simulation of 
biomolecules and/or nano-molecular systems comprising: 

constructing at least one set of equation representing the molecular 
20 dynamics of at least one molecule of the biomolecules and/or the nano- 
molecular systems; 

constructing an electronic circuit representing every set of equation; 

and 

determining molecular dynamics simulation by measuring voltage at 
25 two or more connection points of the circuit. 

The equations may be linear and/or non-linear second order ordinary 
differential equations (ODEs). Preferably, non-linear second order ordinary 
differential equations (ODEs). 
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The biomolecule and/or nano-molecular system of the invention comprise in 
general, but not exclusively, amino acids, nucleotides, organic molecules, 
and/or inorganic molecules. 

According to one aspect of the invention, construction plan is disclosed for 
5 circuit group of each unit in a biomolecule (amino acid or nucleotide) or a 
nano-molecular-system to represent all its internal interaction bonded and 
non-bonded interactions and bonded interactions to its nearest neighbouring 
amino acid or nucleotide. Each circuit group is composed of at least one of: 
(1) internal bond stretch, angle bending, torsion, non-bonded units; and (2) 
10 bond stretch, angle bending, and torsion units to its two nearest neighbouring 
amino acid or nucleotide. In a molecular dynamics simulation of a protein or 
nucleic acid, these circuit groups can be combined together according to the 
sequence of the protein or nucleic acid. Additional non-bonded units 
representing inter-residue non-bonded interactions can be added to each 
15 circuit group accordingly. 

Accordingly, the electronic circuit of the invention may comprise at least one 
atom-position circuit unit, wherein the atom-position circuit unit represents the 
position of an atom of a molecule or a molecular system. The atom-position 

20 circuit unit comprises at least one atom-atom interaction circuit subunit, the 
atom-atom interaction circuit subunit representing a sub-unit of atom-atom 
interactions within a molecule or a molecular system and comprising at least 
one of: internal bond stretch, angle bending, torsion, non-bonded units; bond 
stretch, angle bending, and torsion units; between at least two nearest sub- 

25 unit of a molecule. 

Each atom-atom interaction circuit subunit represents a term in the molecular 
dynamics equation and comprises at least one of the following: bond stretch x 
unit, bond stretch y unit, bond stretch z unit, angle bending x type-A unit, 
angle bending x type-B unit, angle bending y type-A unit, angle bending y 
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type-B unit, angle bending z type-A unit, angle bending z type-B unit, torsion x 
type-A unit, torsion x type-B unit, torsion y type-A unit, torsion y type-B unit, 
torsion z type-A unit, torsion z type-B unit, non-bonded x unit, non-bonded y 
unit, non-bonded z unit, hydrogen-bond x unit, hydrogen-bond y unit, and 
5 hydrogen-bond z unit; and 

wherein x, y, and z represent the coordinates of each atom of the 
molecule, and type-A represents the case of the atom being in the middle- 
position of an angle bending or torsion connection with other atoms, and type- 
B represents the case of the atom being in the end-position of an angle 
10 bending or torsion connection with other atoms. 

According to a particular aspect, all the twenty subunit elements listed above 
may be used for the construction of the circuit units and electronic circuits of 
the invention. 

In a further aspect of this invention, methods are disclosed for maintaining the 
15 voltages in the circuits (xyz coordinates and functions of xyz coordinates) 
within the allowed range. Scaling factors are applied to x, y, z coordinates and 
to the molecular dynamics equations so as to keep the x, y, z coordinates and 
each term in the equations within the allowed range. In addition, resistors and 
amplifiers can be used to scale-down and scale-up the voltages at certain 
20 segment of the circuits. Moreover automatic gain control circuits can also be 
used to ensure the voltages are within the required range. The automatic gain 
control circuits are described in "Function circuits: Design and applications" 
Y.J. Wong and W.E. Ott, McGram-Hill, 1976. 

25 Accordingly, the method further comprises maintaining the voltage level in the 
circuit between two fixed voltage values. In particular, x, y and z represent the 
coordinates of the molecule, and the voltage level of the circuit is maintained 
between two fixed voltage values by: 
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applying (non-zero real numbers) scaling factors to the x, y and z 
coordinates and to the molecular dynamic equation; 

applying at least one resistor and/or amplifier at one or more 
connection point of the circuit, thereby scaling-down or -up the voltage of one 
5 or more segment of the circuit; and/or 

applying automatic gain control circuits. 

According to another aspect of the second embodiment, it is provided a circuit 
group representing the interaction pattern in the chemical structure of a 
10 molecule or a sub-unit of interaction pattern in the chemical structure of a 
molecule comprising at least one of the following: 

a bond stretch connection between each atom pair of the molecule 
covalently bonded to each other; 

an angle bending connection pair between a first atom and other two 

15 atoms; 

a torsion connection bundle between a first atom and other three 
atoms; and 

a non-bonded connection between each atom pair whose atoms are at 
least four bonds away from each other. 

20 

A circuit unit according to the invention may comprise at least one circuit 
group as disclosed. 

The invention also provides an electronic circuit comprising at least one circuit 
25 unit, the circuit unit comprising at least one circuit group, the circuit group 
representing a sub-unit of interaction pattern in the chemical structure of a 
molecule and comprising internal bond stretch, angle bending, torsion, non- 
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bonded units; and/or bond stretch, angle bending, and/or and torsion units; 
between at least two nearest sub-unit of a molecule. 

The electronic circuit represents a term in the molecular dynamic equation, 
5 and wherein the circuit unit comprises at least one of the following: bond 
stretch x unit, bond stretch y unit, bond stretch z unit, angle bending x type-A 
unit, angle bending x type-B unit, angle bending y type-A unit, angle bending y 
type-B unit, angle bending z type-A unit, angle bending z type-B unit, torsion x 
type-A unit, torsion x type-B unit, torsion y type-A unit, torsion y type-B unit, 
10 torsion z type-A unit, torsion z type-B unit, non-bonded x unit, non-bonded y 
unit, non-bonded z unit, hydrogen-bond x unit, hydrogen-bond y unit, and 
hydrogen-bond z unit; and wherein x, y, and z represent the coordinates of 
each atom of the molecule, and type-A unit represents the case of the atom 
being in the middle-position of an angle bending or torsion connection with 
15 other atoms, and type-B represents the case of the atom being in the end- 
position of an angle bending or torsion connection with other atoms. 

According to a particular aspect, the invention relates to design plans for the 
construction of circuit units representing all possible terms in the ODEs for 
molecular simulation of proteins, nucleic acids and organic molecules. The 
20 units are those represented above. 

The invention also provides a method for the manufacture of an electronic 
circuit representing at least one biomolecule and/or nano-molecular system, 
the electronic circuit comprising at least a unit circuit comprising at least a 
circuit group, wherein the circuit group represents a sub-unit of interaction 
25 pattern in the chemical structure of a molecule or a molecular system 
comprising: 

introducing a bond stretch between each atom of a pair of atoms 
covalently bonded to each other; 
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introducing an angle bending connecting pair between a first atom and 
other two atoms; 

introducing a torsion connection bundle between a first atom and other 
three atoms; and 

introducing a non-bonded connection between each atom pair whose 
atoms are at least four bonds away from each other. 



The equations of motion of molecular dynamics simulation of a molecule can 
10 be derived from Newton's second law of motion given the potential energy 
functions, force fields and a starting 3D structure. The potential energy 
functional for proteins, nucleic acids and organic molecules is given by: 

V = 1/2 Ebon* K r (R - Req) 2 + 1/2 Wangles Ke (9 - 6 eq ) 2 +1/2 S, ore ions V„ [ 1 + COS(n<H) ] + 
1 5 E„on bonded [ V^ 2 - + qfli k, ^ + S H bonds V H (0 U ) 

as described in "A second generation force field for simulation of proteins, 
nucleic acids, and organic molecules", W. D. Cornell et. al., J. Am. Chem. 
Soc. 117, 5179-5197 (1995); and in "CHARMM: A program for 
20 macromolecular energy, minimization, and dynamics calculations", B. R. 
Brooks, et.al. J. Comp. Chem. 4, 187-217 (1983). 

The hydrogen bond potential V H (r) is different for different force fields. The 
following is a list of some of the V H (r) appeared in the literatures: 

25 

V H (r) = A/r 12 - B/r 6 + qi^ /s r r 9 for AMBER (same as non-bonded term) 

= (A/r 12 - B/r 10 )cos m (e A .H.D)cos n (e A A.A.H)sw 1 (r)sw 2 (e) for CHARM 
_ v 0 (1-e" a(r " ,0) ) 2 - v o for Prohofsky/Chen (2) 

30 All force field parameters for each atom or atom-atom pair in a biomolecule or 
nano-molecular-machine are given in the literature. For example, AMBER 
force fields are described in "A second generation force field for simulation of 
proteins, nucleic acids, and organic molecules", W. D. Cornell et. al., J. Am. 
Chem. Soc. 117, 5179-5197 (1995). As for organic molecules, all parameters 
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other than partial charges are available. Partial charges of an organic 
molecule can be computed from quantum chemistry software such as 
Gaussian. 

From Newton's second law, the atomic position of the i-th atom, in terms of its 
xyz coordinates (xj, y„ Zj) is given by: 



25 



30 



m i d 2 x i /dt 2 = -dV/dx i 

= -S b0 ndsK r (R-R eq )B r xi - Z angles Ke(e-0 eq )B 9 xi + 2:tor S ionsnV n sin(n<hr) B\, 

- Z nonb onc.ed [6Bi/ rij 7 - IZAi/rij 13 - q iqi /e r r* ](xrXj)/r, - Abends [dV H (r)/dr ](x r x h )/r 

m i d 2 y l /dt 2 = -dV/dy i x . 

= -2bond8Kr(R-Re«0B r y. - ^angles Ke (9 " Oeq) B % + ^torsions nV n Sin(n<j>- Y ) B% 

- Snonbonded [6B^r, 7 - 12A^ 18 " <W|/* r* ]<y r y,yr, - E Hb onds [dV H (r)/dr KyfYhVr 

m i d 2 z i /dt 2 = -dV/dz i i 
= -Stonds K r (R - R eq ) B r 2i - S angle s K, (9 - 6 eq ) B 6 zi +/ S torS ions nV n sin(n<t)-Y) B*» 

- Snonbonded [6B«/r, 7 - 1 ZAi/r.j 13 - qflj/s, n, Kzrzj)/^ - I H bonds [dV H (r)/dr ](z r z h )/r 
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(3) 



(4) 



(5) 



where mi is the mass of the i-th atom; B r xi , B r yi , B r 2i are the x, y z components 
of the B-matrix for bond stretching; B e xi , B e yi . B 9 zi are the x, y, z components 
for the B-matrix for bond angle bending; B* xi , B* yi , B* zi are the x, y, z 
components for the B-matrix for torsion. As described in "Vibrational states", 
S. Califano. John Wiley & Sons, New York (1976), these B-matrix elements 
are given by: 

(6) 
(7) 
(8) 



if I-th atom is in the middle 

if I-th atom is at end (9) 

if I-th atom is in the middle 

if I-th atom is at end (1 0) 



B f x i — 


-Ai] 




B f yi = 


-Bij 




B r zi = 


-Cij 




R 9 - 
O X i - 


Ea = 
Eb 


-(Eb+Ec) 


R e - 

B yj - 


Fa = 


-(Fb+Fc) 



= Fb 
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B 6 zi = Ga = -(Gb+Gc) 
= Gb 



if I-th atom is in the middle 

if I-th atom is at end (1 1 ) 
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B* xi = Ua = 1/rnPa 2 + qcSs/HkPc 2 - qaSe/rikPa 2 if l-th atom is in the middle 

= Ub = - s^/ryPa 2 if l-th atom is at end (12) 

B\, = Va = 1/rijPa 2 + qcS3/r ik Pc 2 - q a sVr ik p a 2 if l-th atom is in the middle 

= Vb = - sl/ryPa 2 if l-th atom is at end (13) 

B* xl = Wa = 1/r iiPa 2 + q c si/r ik p c 2 - qaS 2 /r,kPa 2 if l-th atom is in the middle 

= Wb = - Sa/rijPa 2 if l-th atom is at end (14) 

where ry, Ay, By, C ij( q a , q c , Pa, Pc, Eb, Ec, Fb, Fc, Gb, Gc, Si, s 2 , s 3 , s 4 , s 5 , s 6 
are given by: 

r lj =[(x i -x i ) 2 + (yry i ) 2 +(zrz i ) 2 ] 1 ' 2 < 15 > 

A ij = (x j -x i )/r ij < 16 > 

Bij = (yj-yO/nj (17) 



(19) 
(20) 

1/2 (21) 



q a = cos9 a = AyAjk + BjjBflc + QjCjk 
q c = cosGc = -(AikAki + B^B^ + CuAa) 
p a = sin9 a =(l-qa) 



■ a /-i „ \i« (22) 
p c = sm0 c =(l-q c ) v ' 

Eb = (Aijqa-AaO/rijPa (23) 



Ec = (Aikqa-AiiVrap 



Fc = (B ik qa-By)/r ik pa 



(24) 



Fb = (Byqa-BiO/ryPa (25) 



(26) 



Gb = (Cijqa-CaO/rijPa (27) 

Gc = (Cfcqa-CyVraPa 

5 1 = A ik B kl -B flc A kl 

5 2 = A ik B ij -B ik A ij ( 30 > 

53 = AkjCflc — CkiAjk 

54 = AyCik — CyAjk ^ ^ 
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Ss^BikOd-CflAd ^ 
Sg = BikQj — CyBik (34) 
This set of nonlinear ODEs contain terms of n it (Xi-Xj)*(x k -Xj), (xi-Xj)/^ , 
sin(n<t>-Y) etc., which can be represented by a corresponding set of electronic 
analog circuits composed of known operational amplifiers, function circuits 
and other circuit components. 

Referring now to Figure 10, Figure 11, and Figure 12, there are disclosed 
design plans for construction of electronic circuits representing molecular 
dynamics simulation equations (3), (4) and (5) for each atom in a protein, 
nucleic acid, or organic molecule respectively. For each equation, its circuit 
starts with a point Xj or y or Zj, which is connected to a switch directed to a 
voltage source representing the initial coordinate. This point is also connected 
to other parts of the circuit or other circuits in a combination of one or more of 
the seven cases as illustrated in Figure 10, Figure 11, and Figure 12. In a 
simulation, the time-dependent x or y or z coordinate of the i-th atom can be 
by measurement of voltages at Xi or y or z, . 

Construction Plan for the Specially Designed Circuit Groups Representing All 
the Internal Terms and Bonded terms with nearest neighbors in the ODEs of a 
sub-unit of a biomolecule or in a nano-molecular-machine 

Referring now to Figure 13 to Figure 15, there are disclosed design plans for 
construction of circuit group for each illustrative sub-unit, amino acid, of a 
biomolecule, protein, to represent all its internal terms and bonded terms to its 
nearest neighboring subunits, amino acids. 

Referring now to Figure 17 to Figure 19, there are disclosed design plans for 
construction of circuit group for each sub-unit of an illustrative sub-unit of a 
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nano-molecular-machine in Figure 16 to represent all its internal terms and 
bonded terms to its nearest neighboring subunits. 

In these figures, each box represents the x, y and z circuit sections of each 
atom of an amino acid or nucleotide. An arrow represents a connection 
between two atoms such that the output voltage Xj, yj, or zj from the x, y or z 
circuit section of an atom is directed to a bond-stretch, angle bending or 
torsion unit in the corresponding section of the atom to which the arrow is 
pointed. For clarity purpose, the connection profile for bond stretch, angle 
bending, and torsion is displayed in separate Figures. In reality, these should 
be put together in the same circuit group configuration. Figure 1 3 and Figure 
17 illustrates connections to bond stretch units of each circuit section, Figure 
14 and Figure 18 illustrates connections to angle bending units of each circuit 
section, and Figure 15 and Figure 19 illustrates connections to torsion units of 
each circuit section. 

A circuit group is in general constructed by the following rules: 

From the chemical structure of an amino acid or nucleotide: 

A bond stretch connection is introduced to each atom pair 
covalently bonded to each other; 

An angle bending connection pair are established from atom A to 
two other atoms B and C, in cases that they are linearly bonded as 
A-B-C or B-A-C or B-C-A etc.; 

A torsion connection bundle are constructed from atom A to three 
other atoms B, C and D, in cases that they are linearly bonded as 
A-B-C-D or B-A-C-D or B-C-A-D etc.; 

A non-bonded connection is introduced to each atom pair whose 
atoms are at least 4 bonds away from each other. 
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Construction Plan for the Specially Designed Circuit Units Representing 
Different Terms in an ODE For Molecular Dynamics Simulation 

Referring now to Figure 20 to Figure 40, there are disclosed design plans for 
5 construction of special circuit units representing all possible terms of an 
ordinary differential equation associated with a biological or chemical reaction 
pathway. The operational amplifiers and random noise generator in these 
units are all described in the literature. The inverting amplifier and 
noninverting amplifier are described in "Electronic circuit analysis and design, 
10 2 nd edition" D.A. Neamen, McGraw-Hill, 2001. The multiplier is described in 
"Fundamentals of linear circuits" T.L. Floyd, Macmillan Publishing Co. 1992. 
The exponentiator is described in "Function circuits: Design and applications" 
Y.J. Wong and W.E. Ott, McGram-Hill, 1976. The sine function amplifier is 
described in "Modern operational circuit design" J. I. Smith, Wiley- 
15 Interscience, 1971. The cos' 1 function amplifier can be generated in different 
ways. One example is the use of circuits of polynomials and another example 
is the use of circuits of non-integer exponent as described in "Function 
circuits: Design and applications" Y.J. Wong and W.E. Ott, McGram-Hill, 
1976. 

20 

Some circuit units in Figure 10 to Figure 40 use one or more of the following 
specially introduced circuit sun-units: rABC unit, AqA unit, pq unit, Ua unit, ss 
unit, si unit, <J> unit, sincj> unit. These sub-units are given in Figure 41 to Figure 
44. 

25 

Methods for Maintaining the Voltages in the Circuits Within Allowed Range 

In large proteins, the distance between two atoms can be larger than a few 
30 hundred angstroms. Thus the range of the xyz coordinates in these proteins 
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may exceed the allowed voltages in a circuit. In addition, some terms in the 
molecular dynamics simulation equations can become very large if two atoms 
get too close or too far apart from each other, thereby causing the voltages at 
certain parts of the circuits to exceed the allowed range. 

These problems can be solved by maintaining the voltage level in the circuit 
between two fixed voltage values, by using at least one of the following 
methods: 

(1) Applying scaling factors to x, y, z coordinates and to the molecular 
dynamics equations so as to keep the x, y, z coordinates and each 
term in the equations within the allowed range. 

This can be accomplished by multiplying one or two appropriately 
selected small numbers (non-zero real numbers) 5 or 6e to the 
Equation 3, Equation 4, and Equation 5. 



(2) Resistors are used to scale-down the voltages at the input end of a 
segment of the circuits and amplifiers are used at the output end of the 
segment to scale-up the voltages to the proper level. 

(3) Automatic gain control circuits can be added to the circuits to ensure 
the voltages are within the required range. 



25 

Illustrative example: Electronic protein emulator 

Structural optimization and molecular dynamics of a protein can be described 
by a set of ODEs. An illustrative amino acid is governed by the collection of 
30 ODEs of each constituent atom (Equation 3-5). The circuit for each atom is 
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illustrated in Figure 10 to Figure 12. The circuit group for all the atoms of an 
amino acid is illustrated in Figure 13 to Figure 15. 

Illustrative example: Electronic nano-molecular-system emulator 

There have been efforts for developing nano-molecular-systems. One 

example is a mimic of an enzyme ribonuclease shown in Figure 16. 

Structural optimization and molecular dynamics of such a nano-molecular- 
machine can be described by a set of ODEs. An illustrative sub-unit of 
cyclodextrin shown in Figure 16(a) is governed by the collection of ODEs of 
each constituent atom (Equation 3-5). The circuit for each atom is illustrated in 
Figure 10 to Figure 12. The circuit group for all the atoms of a sub-unit is 
illustrated in Figure 17 to Figure 19. 



